Computing stationary distributions in equilibrium and 
non-equilibrium systems with Forward Flux Sampling 



Chantal Valeriani, 1 Rosalind J. Allen, 2 ' 1 Marco J. 
Morelli, 3 Daan Frenkel, 3 and Pieter Rein ten Wolde 3 

1 FOM Institute for Atomic and Molecular Physics, 
Kruislaan 407, 1098 SJ Amsterdam, The Netherlands. 
2 SUPA, School of Physics, The University of Edinburgh, 
Mayfield Road, Edinburgh EH9 3JZ, UK. 
3 FOM Institute for Atomic and Molecular Physics, 
Kruislaan 407, 1098 SJ Amsterdam, The Netherlands. 

We present a method for computing stationary distributions for activated processes 
in equilibrium and non-equilibrium systems using Forward Flux Sampling (FFS). In 
this method, the stationary distributions are obtained directly from the rate constant 
calculations for the forward and backward reactions; there is no need to perform 
separate calculations for the stationary distribution and the rate constant. We apply 
the method to the non-equilibrium rare event problem proposed by Maier and Stein, 
to nucleation in a 2-dimensional Ising system, and to the flipping of a genetic switch. 

I. INTRODUCTION 

Rare events are ubiquitous in physics, chemistry, and biology; examples include crystal 
nucleation, chemical reactions, and protein folding. Rare events are activated processes, 
for which the average waiting time between events can be orders of magnitude longer than 
the duration of the event itself. This makes these events intrinsically difficult to investigate 
experimentally. Computer simulations are therefore a natural tool to use - yet conventional 
numerical techniques are impractical for rare events, because most of the CPU time is wasted 
on the uneventful waiting time between events. A number of "rare event" simulation schemes 
have recently been developed in the field of soft-condensed matter physics, which make it 
possible to zoom in on the rare events themselves. Techniques such as umbrella sampling 
allow the calculation of free-energy barriers separating the stable states , while 



schemes such the Bennet-Chandler method jf], [f| also allow the computation of rate con- 
stants. Transition path sampling 0, Q, S, Q| allows both rate constants and transition paths 
to be obtained. These techniques have been used for a wide range of applications including 
ion permeation through membranes, protein folding, and nucleation. However, these schemes 
require prior knowledge of the phase-space density. For systems that are in thermodynamic 
equilibrium — with detailed balance and microscopic reversibility — the phase-space density 
is known: it is given by the Boltzmann distribution. In contrast, for systems that are out of 
equilibrium, the phase-space density is usually not known. This means that most numerical 
techniques for simulating rare events are limited to equilibrium systems, and thus exclude 
a host of important rare-event problems in non-equilibrium systems, such as polymer col- 
lapse under flow, crystal nucleation under shear, and rare events in biology, such as protein 
translocation and switching events in biochemical networks. We have recently developed a 
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numerical technique, called Forward Flux Sampling (FFS)[lll. 1 12. Il3j. that makes it possible 
to compute rate constants in both equilibrium and non-equilibrium systems with stochastic 
dynamics. In this paper, we show how stationary distributions can also be obtained directly 
from an FFS calculation, for both equilibrium and non-equilibrium systems. For equilib- 
rium systems the advantage is that from an FFS simulation one can obtain not only the 
rate constant, but also information about the free-energy landscape. For non-equilibrium 
systems the concept of free energy does not apply, but one can obtain the steady-state prob- 
ability distribution as a function of a chosen order parameter (or order parameters). To our 
knowledge, this is the first method to be proposed for computing stationary distributions 
for multi-dimensional non-equilibrium systems that are in steady state. 

In soft-condensed matter physics, the rate k of an activated process in an equilibrium 
system is often written as the product of two factors: 

k = R(q*)p(q*). (1) 

Here, q is an order parameter that connects the initial and final states, assuming that the 
system evolves between two states. It is defined such that for q < q* , the system is in the 
initial state, while for q > q* it is in the final state. The quantity p(q*) is the probability 
that the system is at the dividing surface q = q* and R(q*) is the rate at which this dividing 
surface is crossed. For equilibrium systems, p(q*) is proportional to exp(— (3AG(q*)), where 
P is the inverse temperature and AG(q) is the (Landau) free energy of the system as a 



function of the order parameter q. It is natural to locate the dividing surface q* at the top 
of the free-energy barrier AG(q) separating the two states. The rate constant is thus given 
by the probability of being at the top of the free-energy barrier, multiplied by a kinetic 
prefactor. 

The Bennett-Chandler method for computing rate constants for activated processes uses 
a two-step procedure jf], one first computes the free-energy barrier, using, for example, 
the umbrella sampling scheme , and then the kinetic prefactor, using a Molecular 

Dynamics simulation in which trajectories are fired from the top of the free-energy barrier. 
However, this method is computationally demanding, and its success depends strongly on 
the choice of the reaction co-ordinate q. If q is poorly chosen, the system will sample the 
wrong part of the phase space, which will not only conceal the mechanism of the transition, 
but also impede the computation of the rate constant — while the choice of q co-ordinate 
does not affect the value of the rate constant k, it can strongly affect the efficiency with 
which k is computed. For high-dimensional complex systems it can be difficult to make a 
good choice for q, since this requires a priori insight into the reaction mechanism. 
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Transition-path sampling (TPS) has been developed to alleviate these problems [7), [8|, LLJ, 



151 ] . This scheme generates an ensemble of trajectories between the initial and final states 
using Monte Carlo sampling in trajectory space. TPS only requires an order parameter to 
distinguish the initial and final states; this order parameter does not need to be the true 
reaction co-ordinate. TPS thus makes it possible to compute the rate constant without 
prior knowledge of the reaction mechanism. However, this method does require knowledge 
of the steady-state phase space distribution, which is needed for the acceptance/rejection 
step in the Monte-Carlo scheme, and it does not allow direct computation of the free-energy 
barrier separating the two states. Moroni et al. have developed a related method, transition 
interface sampling (TIS), which relies on the computation of crossing probabilities of a series 



of interfaces between the initial and final states (si, Q, lid, [l7|. A variant of this method, 
partial path TIS (PPTIS), assumes loss of time correlations in the transition paths over a 
distance of two interfaces. Moroni et al. have recently shown how the free-energy barrier as 
well as the rate constant can be obtained from a single TIS / PPTIS calculation fl8|. As 
for TPS, both TIS and PPTIS require the system to be in thermodynamic equilibrium. The 
"milestoning" method of Faradjian and Elber also employs a series of interfaces to compute 
rate constants and also assumes that the interface-crossing probability does not depend 



upon the full history of the path {19I ]. In an alternative approach, Vanden-Eijnden et al have 
developed a set of "string" methods, which can be used to compute minimum free-energy 



paths and the probability current of reactive trajectories for equilibrium systems [20l. I2l|. 

The algorithms discussed above — TPS, (PP)TIS, milestoning and the string methods — 
are limited to systems that are in thermodynamic equilibrium. The Forward Flux Sampling 
(FFS) method, and its variants, were developed to calculate rate constants and transition 
jaths for rare events in equilibrium and non-equilibrium systems with stochastic dynamics 
111, [l2, Like TIS, PPTIS, and milestoning, FFS uses a series of interfaces to compute 
the rate constant. However (unlike PPTIS and milestoning), FFS does not make the Marko- 
vian assumption that the distribution of paths at the interfaces is independent of the path 
histories. The order parameter that is used to define the location of the interfaces need not 
be the reaction co-ordinate, and the choice of order parameter, in principle, does not bias 
the dynamics of the transition paths. 

We have recently shown that the rate constant for activated processes in non-equilibrium 



systems that are in steady state can also be written in the form of Eq. [T] [22 ] . The quantity 
p(q) is then the stationary probability distribution function for the order parameter q. In this 
paper we show that the stationary distribution p(q), as well as the forward and backward rate 
constants and transition paths, can be obtained by performing two FFS calculations — one 
for the transition from the initial to the final state, and the other for the reverse transition. 
The method can be applied to both non-equilibrium and equilibrium systems; in the latter 
case, p(q) corresponds to the Boltzmann distribution. The method is conceptually similar 
to that used in TIS and PPTIS to compute free-energy barriers, in the sense that the 
stationary distribution p(q) is obtained by matching the forward and backward trajectories 
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In the next section, we explain the FFS algorithm . In sections [TTT] and HVl we 
discuss the theory and method for obtaining stationary distributions. We then illustrate the 
method using symmetric and asymmetric double- well potentials (section IVTl) . and the two- 
dimensional non-equilibrium rare event problem proposed by Maier and Stein (section lYlIj ). 
In section IVHH we use the method to calculate the free-energy barrier for nucleation in a 
two dimensional Ising system. Finally, in section ILXl we compute non-equilibrium stationary 
probability distributions for a bistable model genetic switch. 



II. FORWARD FLUX SAMPLING 



We consider rare, spontaneous transitions between two regions of state space A and B. 
The phase space co-ordinates are denoted by x and the regions A and B are defined in terms 
of an order parameter X(x) such that the system is in state A if A(x) < Ao, and it is in state 
B if A(x) > A n . The key principle is to use a series of interfaces A , Ai, . . . , A n _i, A n , to drive 
the system from state A to state B in a ratchet-like manner. The idea of the interfaces is 
that they make it possible to capitalize on all those fluctuations that bring the system in 
the direction of the final state B. 

Supposing that with a conventional (say MD) simulation, the system exhibits a rare 
fluctuation that moves it up the barrier, and that it crosses an interface between state A 
and the top of the barrier, if we would continue this succesfull run, then most likely the 
system would roll back down the hill, i.e. relax back towards state A, and one would have to 
wait for "another" rare fluctuation that moves the system in the direction of B. By storing 
the configurations at the interfaces, we can thus efficiently exploit all those fluctuations that 
move the system up the barrier. 



In this paper, we make use of the original FFS scheme 
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Ill ] presented in more details in 
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In FFS, one first performs a conventional, brute-force simulation in state A. Each time 
the system crosses the interface Ao in the direction of increasing A during this simulation, 
the co-ordinates of that state point are stored. One also measures the average number per 
unit time of these crossings. At the end of this simulation, one has a measure of the flux 
$a of trajectories crossing Ao from A, as well as a collection of state points corresponding 
to crossings of the first interface, A , coming from A. This collection is then used to provide 
starting points for a set of trajectories, each of which is continued until the system either 
reaches the next interface, Ai, or returns to state A (i.e. reaches A ). This procedure 
generates a new collection of state points at the next interface, which are the end points of 
those trajectories that arrived at Ai from Ao- One also obtains an estimate of the probability 
-P(Ai|Ao) that a trajectory which reaches Ao from A will subsequently reach Ai without 
returning to A - this is simply the fraction of trajectories which arrive at Ai. By repeating 
this procedure for all subsequent interfaces, one has for each interface i an estimate of the 
probability P(Aj + i|Aj) that, given that a trajectory has reached interface i coming from A, 



it subsequently reaches A i+ i before returning to A. The rate constant k A B can then be 
obtained from fis| 

n— 1 

k AB = $ A l[P(\ l+1 \\ i ). (2) 

i=0 

By tracing back paths that successfully arrive at A n , one can also sample the transition path 
ensemble for the rare event. Analysis of these paths can lead to insight into the mechanism 
by which the event occurs. 



III. STATIONARY DISTRIBUTIONS: THEORY 

We are interested in computing the stationary distribution p(q), where p(q)dq is the 
probability of observing the order parameter q in the range q — > q + dq, for a system that is 
in a stationary state. We stress the fact that the order parameter q for the computation of 
the stationary distribution function need not be the same as the order parameter A that is 
chosen for the FFS calculation. The stationary distribution can be expressed as 

p(q) = (S(q-q(x))), (3) 

where x is a point in the multi-dimensional phase space. For equilibrium systems, the 
contributions to the average in Eq. [3] are weighted according to the Boltzmann distribution, 
while for non-equilibrium systems they are weighted according to the steady-state phase- 
space density. For both equilibrium and non-equilibrium systems that are in steady state 
and ergodic, this ensemble average is equivalent to a time average over a long brute-force 
simulation, in which p(q) measures the frequency with which value q of the order parameter 
is "visited" by the trajectory. 

The distribution function p(q) is easy to sample close to the stable states A and B, using 
conventional, brute-force simulation. However, this method will lead to poor statistics in 
the "barrier" region between A and B, which is rarely visited. We use FFS to obtain p(q) 
in the "barrier" region, and supplement this with conventional sampling in the two stable 
states to obtain the complete distribution function. 

The key idea which we use to obtain stationary distributions with FFS is to divide the 
"visits" of an imaginary, very long simulation trajectory to value q of the order parameter 
into two categories, according to whether the trajectory was most recently in state A or 



state B. We therefore write p(q) as the sum of two contributions 

p{q) = ^A{q) + ^ B {q), (4) 

where i/)a(q) is the contribution to the probability density p(q) from those trajectories that 
come from region A, and ipsio) is the contribution due to trajectories coming from B (see 
Fig. ED. 
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Figure 1: A sketch of all the possible trajectories that contribute to the stationary distribution p(q) 
(see Eq. [3]): Trajectory 1 comes from A and goes back to A; trajectory 2 comes from A and goes 
to B; trajectory 3 comes from B and goes back to B; trajectory 4 comes from B and goes to A. 
The FFS simulation from A to B harvests the trajectories corresponding to types 1 and 2, while 
the FFS simulation for the reverse transition generates trajectories of types 3 and 4. The interfaces 
{Ao, Ai, . . . , A n _i, A n } used in the FFS simulation are also shown. 

In the basins of attraction A and B, the trajectories will quickly lose memory of where 
they came from - i.e. we expect excursions out of a basin of attraction to be uncorrelated. 
This, as we describe below, makes it possible to obtain the distribution function Vu(?) from 
an FFS simulation for the transition from A to B, while ipsio) can be computed using an 
FFS simulation for the reverse transition (see Fig. []]). For equilibrium systems, the free 
energy profile can be obtained from AG ~ — k B T\n [p(q)] once p(q) is known. 

The function i/)a(q) is given by 

Mq) = PA®AT+(q;\ ). (5) 

Here, pa is the probability that the system is in state A and $^ is, as in Eq. [2, the flux 
of trajectories leaving state A (i.e. crossing the surface A coming from A). The quantity 



r + (g; Ao) = (S(q — q(x)))\ is the average time spent at order parameter q by a trajectory 
that originates from interface A . We note that r + (g; A ) includes contributions both from 
paths that start in A and ultimately reach B, and from those that start in A and ultimately 
return to A without reaching B (see Fig. CD). 

In FFS, we use a series of interfaces to sample the phase space between A and B in stages. 
At each stage, an ensemble of paths is generated by firing off trajectories from points on 
an interface that have been obtained in the previous stage; each of these trajectories is 
terminated as soon as it reaches either the next interface or A (see section HI!)- We denote 
the average time spent at order parameter q, for a trial run that is fired from interface Aj (and 
terminated at Aj+i or Ao) in the FFS procedure, by n + (q;Xi). As shown in the appendix, 
r + (q; Ao) is then given by: 

n— 1 i— 1 

T+(q;X ) = vr+(g;A ) + ^7r+(g; Xi) JJP(A i+ i|A i ). (6) 

i=l j=0 

The factor [TV=o P ( A i+il A i) 

reweights the distribution 7r(g; A«) to correct for the enhanced 
sampling at interface i which has been achieved by the FFS procedure. This factor is a 
direct output from the FFS simulation (see section [III Eq. El and the appendix). The FFS 
calculation for the forward transition thus yields k^B, and r + (g; A ). 

To calculate p(q), we also need to evaluate i[>B(q) in Eq. 31 by carrying out an FFS 
calculation in the reverse direction, from B to A. The entire FFS algorithm is carried out 
in reverse: in the initial, brute-force simulation, we begin with the system in state B and 
collect crossings of interface X n coming from B. We fire trajectories from Xi which either 
reach Aj_i or return to X n . The result is a value for the reverse rate constant ksA, the 
flux and the distribution functions 7r_(g; Xi) for the order parameter q, sampled over the 
ensemble of paths that are fired from interface A, and terminated at Aj_i or A n in the reverse 
FFS procedure. These are related to the distribution function r_(g; A„) for all trajectories 
leaving A n from B by: 

1 i+l 

r_(g; A n ) = 7r_(g; X n ) + vr_(g; A;) JJ P(A i _ 1 |A J ), (7) 

i=n—l j=n 

where P(Aj_i|Aj) are the conditional probabilities of reaching interface j — 1 from Xj, eval- 
uated in the reverse FFS procedure. The distribution i/jb(q) is then given by: 

Mq) = PB®BT-{q;\n). (8) 



To obtain pa and ps in Eqs([5|) and (JSj) , we note that in steady state 



Pa^ab = Pb^bAi (9) 

where Icab and kg a are the forward and backward rate constants measured in the forward 
and backward FFS calculations, respectively. Since we are assuming a two state system (i.e. 
ignoring intermediate states), we also know that Pa+Pb = 1. This implies that 

k B A/k A B , n v 

1 + ksA/kAB 

and 

1 + k B A/k A B 

Combining all this information and using EqpJ, we can obtain the stationary distribution 
function p(q) in the region Ao < A < A n . This can be combined with brute-force sampling 
in the A and B basins to determine p(q) over the full range of q values, if required. 



IV. STATIONARY DISTRIBUTIONS: METHOD 

As discussed above, to obtain the stationary distribution p(q) in the region A < A < A n 
we perform one FFS simulation for the transition from A to B and one for the reverse 
transition. For details on the implementation of the FFS method to com put e the fluxes 
$,4 and $b, as well as the rate constants kAB and ksA, we refer to ref. Hfl. Here, we 
briefly discuss how r + (q; Ao) and r_(g; A n ) are obtained in practice. We consider r + (q; Ao); 
r_(g; A n ) is obtained similarly, but in reverse, as described above. Our aim is to calculate 
the quantities vr + (g; Aj) and P(Xj + i\Xj) in Eq. [6] [or alternatively for the reverse transition, 
7r_(g;Aj) and P(Aj_i|Aj) in Eq. [7]. Considering only the forward FFS procedure: at each 
interface Aj we fire a total of Mj trial runs, each of which is terminated when the system 
reaches either Aj + i or Ao- The probability P(Aj + i|Aj) is then estimated as 

^ilAi)^, (12) 

where Nf is the number of trials that have successfully reached Aj+i. The function n + (q; A«) 
is given by 

iV 

*+<^> = A^? (13 » 



where N q is the number of times that during this set of trial runs the order parameter of 
the system has a value between q and q + Aq. This is given by N q = At J2k=o Ss=o hq( x k,s)i 
where the double sum runs over all the n k steps of all the Mj trial paths starting at interface 
Aj and h q (x) is an indicator function that is one if during a time step the system is between 
q and q + Aq, and zero otherwise; again, note that n k varies from one path to the next. 
The simulation timestep At can in fact be neglected, since is it a constant and we plan to 
normalise p(q) in any case. For algorithms in which the time step can vary, N q is given 
by N q = Yl!k=oYTs=o^tk, s h q (xk,s)-: where At kji is the magnitude of time step s of path k. 
To obtain r + (q; X ) we reweight TT + (q; A«) and sum over all interfaces using Eq. El Once 
i^a(q) and ipB{q) have been obtained by performing FFS simulations in both directions, 
p(q) can be obtained via Eq. [3j We note that ^a{q) and ^b(q) should not be individually 
normalized, since they are not probability distribution functions in their own right, but 
simply contributions to the distribution function p(q). If the average path length for paths 
originating in A and B is different, then the integral of ij)A{q) and ipB(q) over q will be 
different. Normalizing ipA{q) and ipB{q) will result in incorrect relative contributions to p(q) 
from trajectories originating in A and in B. We also note that, when evaluating N q , it is 
important not to double-count the start and end points of trial runs - if the initial point of 
a trial run is deemed to count towards the N q histogram for that interface, then the final 
point should not count as it will be counted as an initial point in the histogram for the next 
interface. 

The above procedure generates p(q) in the region Ao < A < A n . To obtain the full 
distribution p(q), we sample using conventional, brute-force simulation the steady-state 
distribution for the order parameter q in the A and B regions. This will result in distributions 
for A < A + AA (A region) and A > A n — AA (B region), where AA is a small overlap. An 
easy way to fit these curves together is to take their logarithms: the three overlapping parts 
for log p(q) can then be fitted together by a least squares fitting procedure (since a constant 
may be added to each without affecting the distribution). The resulting full profile p(q) is 
obtained by exponentiating log p(q) , and the stationary probability distribution can finally 
be normalised. 



V. STATIONARY DISTRIBUTIONS OF MULTIPLE ORDER PARAMETERS 



It is important to point out that the procedure described in section [TV] may be adapted 
to allow the computation of stationary distribution functions of several order parameters 
("free energy landscapes" in the equilibrium case). In the case where we wish to find the 
stationary distribution (for A < A < A n ) as a function of two order parameters q and r, 
Eq. [H is replaced by 

p(q,r) =ip A (q,r) +t/; B (q,r), (14) 

where 

Mq,r) = p A ^AT + (q,r; A ) (15) 
^B{q,r) = PB$BT-(q,r;\ n ) 

and 

n— 1 i— 1 

T+(q,r;\ ) = vr + (g,r; A ) + ^vr + (g,r; Ai) JJP(A i+ i|A j ) (16) 

8=1 j=0 
1 i+l 

T-(q,r;X n ) = vr„(g, r; A„) + vr_(g, r; A») JJ P(A j _i|A i ). 

i=n— 1 i=" 

To evaluate the functions 7r + (g, r; Aj) and 7r_(g, r; Aj), we use a two-dimensional histogram 
N qr in the co-ordinates q and r: 

AtiV or , x 

= AqArWi < 17) 

and the equivalent for the reverse FFS procedure. Here, N qr is the number of timesteps 
during the set of Mi trial runs fired from Aj for which the system has a value of q between 
q and q + Aq and a value of r between r and r + Ar. 

VI. TESTING ON A ONE-DIMENSIONAL SYSTEM 

As an initial test, we have applied the method to a single particle moving with Brownian 
dynamics in a one-dimensional double-well potential 

V(x) = -bx 2 + ex 4 , (18) 



with b = 2 and c = 1. Distance is measured here in units of Xq, while time is measured in 
units of to. The stationary distribution function, as a function of the x-co-ordinate, is the 
Boltzmann distribution: 

p(x) ~ e- v ^' kBT . (19) 

The system is symmetric, so that that p(A) = p(B). The particle moves according to: 

D 



v(t) 



f(t)+m, 



(20) 



where / is the instantaneous force, D is the diffusion constant and £ is chosen at random 



from a Gaussian distribution with zero mean and variance (£ 2 ) = 2Ddt 231]. We use the 
following values: D = O-Olxg/to, k B T = 0.1 and dt = 0.05to- We have carried out FFS 
simulations with n = 8 interfaces, Ni = 10000 points at interface A , and parameters as 
shown in Table fl] 
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Table I: Interfaces and the number of trials at each interface for the FFS sampling of the symmetric 
one dimensional double-well potential. 



We obtained a forward rate constant Uab — 3.87 ± 0.05 x 10~ 6 tQ 1 (repeating twice to 
obtain error bar). Because of the symmetry of the problem, it was not necessary to carry out 
separate FFS calculations for the forward and backward transitions in this case - the back- 
ward probability distribution can be obtained from the forward one by a simple co-ordinate 
inversion. The stationary distribution obtained from the FFS calculation is compared to the 
expected Boltzmann distribution in Fig. [2l 

We have also considered the asymmetric case, in which a term linear in x is included in 
Eq.UHl 

V{x) = ax + -bx 2 + cx 4 , (21) 

with a = 0.25, b = 2, c = 1, D = 0.01xg/* , k B T = 0.1 and dt = 0.05t . In this case, 
p(A) 7^ p(B) and it is necessary to carry out FFS sampling in both directions. We carry out 



100 




Figure 2: Stationary distribution (solid line) obtained using the procedure described above, com- 
pared to the normalized Boltzmann distribution (circles) for a symmetric double well potential. 
The dotted and dashed lines show iPa(x) and iPb(x) respectively. 

FFS simulations, again with n = 8 interfaces and N% = 10000. For the forward transition, we 
used A = x, and for the backward transition, A = — x. For both the forward and backward 
transitions, the parameters for the FFS runs were as shown in Table HT1 
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Table II: Interfaces and the number of trials at each interface for the FFS sampling of the asymmetric 
one dimensional double-well potential. 

The forward and backward rate constants were calculated to be Uab = 3.03 ± 0.06 x 
lO" 7 ^ 1 and k BA = 3.96 ± 0.03 x lO -5 ^ 1 , and the fluxes across the A boundary were 
$4 = 0.1526 ± 0.0007^ and $ B = 0.3648 ± O.OOOltQ 1 , respectively. Fig. Eh shows r+(x; A ) 
and r_(a;;A n ), while Fig. [3b shows p(x), calculated from Eq. [3] and normalized. Excellent 
agreement is obtained with the expected Boltzmann distribution. 
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Figure 3: (a) Computed stationary distribution p(x, y) as a function of x for y = —0.39 (solid line), 
y = —0.19 (dashed line) and y = 0.01 (dot-dashed line), compared to the expected Boltzmann 
distribution (indicated by circles) for the Maier-Stein system with e = 0.1 and a = /J, = 1. (b) 
p(x,y) as a function of y for x = —0.312 (solid line), x = —0.152 (dashed line) and x = 0.008 
(dot-dashed line). 



VII. TESTING ON THE TWO-DIMENSIONAL MAIER-STEIN SYSTEM 



We now move on demonstrate the calculation of two-dimensional stationary distribu- 
tions using a rare event problem in two dimensions that may be in or out of equilibrium 



overdamped Brownian motion in the force field proposed by Maier and Stein 
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26|]: 



x = f x (x,t) +£ x (t) 
V = /vCM) + f»(t), 



(22) 



where x = (x,y). The force field f = (f x , f y ) (which is time-independent) is given by: 

f x = x — x 3 — axy 2 (23) 
fy = -/^/(l + £ 2 ) 

and the stochastic force £ = (£ x ,£y) results from 5-function-correlated white noise with 
variance e, such that 

(&(*)> = 0; W)£ j (0)) = e6 ij ,6t, (24) 

where i — x,y. This system is bistable, with stable points at (±1,0) and a saddle point at 
(0,0). When a = //, the force field can be expressed as the gradient of a potential energy 
function and the system can be considered to be "at equilibrium", while when a^/J, the force 
field f cannot be expressed as the gradient of a potential and the system is thus intrinsically 
non-equilibrium. In these simulations, we use e = 0.1. Taking A = x, we follow the 
procedure described in section |V] to calculate the stationary distribution for —0.8 < x < 0.8 
as a function of the two order parameters x and y. For the FFS calculations, we use 8 
interfaces, A = —0.8 and A 7 = 0.8, and Ni = 100000 initial configurations at A . The 
parameters used are listed in Table IIII1 
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Xi 


M{\i) 


i 


A? 


M(\i) 





-0. 


1000000 


4 


0.0 


200000 


1 


-0.6 


500000 


5 


0.2 


120000 


2 


-0.4 


300000 


6 


0.4 


100000 


3 


-0.2 


250000 


7 


0.6 


100000 



Table III: Interfaces and the number of trials per interface for the Maier-Stein system. 

We initially consider an equilibrium case, with a — /J, — 1. In this case, the particle moves 
in the potential field 4>(x,y) = y ^ x ^ — ^ + T - Figures S](a) and (b) show the stationary 
distribution p(x,y), as a function of x for y = —0.39, y = —0.19 and y = 0.01, and as a 
function of y for x = —0.312, x = —0.152 and x = 0.008. In both panels, the results are in 
excellent agreement with the expected Boltzmann distribution (shown by circles). 

We next discuss the non-equilibrium case (a ^ //), taking a = 6.67, \i = 2.0 and e = 
0.1. Fig. [5] shows equivalent results to Fig. [H but this time the FFS results are compared 
to stationary distributions computed from long brute-force simulations. The brute-force 




(b) y 

Figure 4: Asymmetric double well potential (a): t + (x; Ao)/Ai (dotted line) and r_(x;A n )/At 
(dashed line) (b): Final result for p{x) obtained from Eq. [3] (solid line) compared to the expected 
Boltzmann distribution (circles). 

simulation results are normalised over all space; the FFS results are multiplied by a constant 
scaling factor to bring them into agreement since they are a priori normalised over the region 
—0.8 < x < 0.8 only. Very good agreement is observed. 

VIII. HOMOGENEOUS NUCLEATION IN A TWO DIMENSIONAL ISING 

MODEL 

We now address a rare event problem in a more complex system: homogeneous nucleation 
in a two-dimensional Ising model. For now, we confine ourselves to an equilibrium system 
without any external shear; non-equilibrium nucleation in an Ising model with an external 




0.01 



0.001 



(b) 




Figure 5: p(x,y) for the Maier-Stein system with a = 6.67, p = 2 and e = 0.1. (a): p(x,y) as a 
function of x for y = —0.39 (solid line), y = —0.19 (dashed line) and y = 0.01 (dot-dashed line), 
(b): p(x,y) as a function of y for x = —0.312 (solid line), x = —0.152 (dashed line) and x = 0.008 
(dot-dashed line). The results of long brute force simulations are indicated by circles. 



shering field will be considered in future work 



27|] . The two-dimensional Ising model consists 



of an L x L square lattice of spins with nearest neighbour interactions and periodic boundary 



conditions. Its Hamiltonian 



28fl 



H 



-J J2 ~ h J2 



(25) 



where J is the coupling constant between neighboring spins {pi = ±1) and h the external 
magnetic field. The prime indicates a sum over first nearest neighbour interactions only. We 
simulate a system with iV = 45 x 45 = 2025 spins, a positive magnetic field (3h = 0.05 and a 
positive coupling constant (3 J = 0.65, above the critical coupling J c . The thermodynamically 



stable state is therefore a ferromagnetic one with net positive magnetization, meaning that 
the system tends to have the majority of its spins in the "up state". However, the state 
with an overall negative magnetization (i.e. spins predominantly in the down state) is 
metastable and the system will remain in that state for a significant time if initialised with 
predominantly down spins. We aim to compute the free-energy barrier, as well as the rate 
constant, for transitions from the metastable "down state" to the thermodynamically stable 
"up state". We begin our simulations in the "down state" and consider the formation of 
a cluster of up spins, under conditions of moderate supersaturation (these conditions are 
identical to those used by Sear [29]). All of our simulations are performed using a Metropolis 
Monte Carlo algorithm, in which we attempt to flip each spin once, on average, during each 
Monte Carlo cycle. 



According to Classical Nucleation Theory 



301 ] . the free energy cost of forming a square 



nucleus of edge length L is given by the sum of a line energy and a surface energy: 

AG = 4 7 L - 2/iL 2 , (26) 

where 7 is the interfacial free energy, h is the driving force for nucleation (magnetic field), 
and —2hL 2 is the energy cost of flipping the whole square nucleus with area L 2 . Using 
Eq. EH the nucleation free energy barrier height is given by 

2 7 2 

AG* = -j-. (27) 



29j, |aij, the 



Plugging in numbers, if we take the interfacial free energy to be /3 7 = 0.74 
barrier height as predicted by classical nucleation theory is [3/S.G* ~ 22. 

We have computed the nucleation free energy barrier using two simulation techniques: 

nnnri 

umbrella sampling [U, |2j, |3J, U| and FFS. In both cases, we characterize the extent of the 
transition using a global order parameter, q = S, the total number of up spins in the 
system. The free-energy barrier is then defined as (3AG(S) = — \a\p(S)/N], where p(S) is 
the probability of observing S up spins in the stationary state. 



For our umbrella sampling calculations, we use a series o 



"windows", defined by a har- 



(JAB,! 



We use 25 windows 



monic potential in S, to bias the sampling of phase space 
to cover the range < S < 300, with an overlap of 11 between neighbouring windows. We 
sample each window for 500000 MC cycles, and fit the resulting histograms together using 
a least-squares fitting procedure to obtain the free-energy profile in the range < S < 300. 



We do not attempt to calculate the barrier for values of S greater than 300, since once the 
top of the barrier is crossed, the system is expected to evolve rapidly and we cannot reply 
on the assumption of local thermodynamic equilibrium. Moreoever, when S is large, the 
growing nucleus is likely to interact with its periodic images in neighbouring cells, making 
the results highly system-size dependent. 

The interfaces Aj for the FFS calculations are also defined in terms of the order parameter 
S. To calculate the free-energy barrier using FFS, we need to be able to sample the reverse 
transition, from the thermodynamically stable "up state" to the metastable "down state". 
In general, this is very difficult for a nucleation problem, since the thermodynamic state is 
much more stable than the metastable state and there is a very high free-energy barrier for 
the system to return to the "down state", making the reverse transition difficult to sample, 
even with FFS. We have overcome this problem in this case by constructing a reflecting wall 
beyond the top of the nucleation barrier. This wall is incorporated via a constraint on the 
system dynamics: each trial move that leads to S > Sb> is simply rejected. Since we are 
only interested in the free-energy profile in the region between A and the top of the barrier, 
we may perturb the free-energy landscape outside this region as we choose. This fact, which 
is also exploited in umbrella sampling, depends on the system being in equilibrium - for a 
driven system, we would not be able to use this approach. The reflecting wall, located at 
S = S B > = 1050, replaces the B state by an artificial stable state B' (see Fig. [6]). 




Figure 6: A schematic view of the free energy landscape. A is the metastable "down state", B is 
the "real" thermodynamically stable state, and B' is the "artificial" stable state, constructed by 
introducing a reflecting wall at S = S B > = 1050. 



The free-energy barrier for the B' — > A transition is much lower than that for the B — > A 
transition, but the shape of the free-energy barrier on the A side remains unchanged. The 
use of the reflecting wall greatly facilitates the FFS calculation for the reverse transition — it 
is possible to carry out the reverse FFS calculations without the wall, but this is rather 
laborious as it requires a large number of interfaces. We have verified that the location of 
the reflecting wall is indeed well beyond the top of the free-energy barrier, which is estimated 
to be at 200 < S < 280. 

State A is defined by the first interface A = 30 - i.e. when < S < 30 the system is in 
the A state. State B' is defined by X n = 1000 - i.e. when 1000 < S < 1050 the system is 
in the B' state. For our FFS calculations, we consider N% = 50 configurations at the first 
interface. The interfaces are located, both for the forward and backward sampling, at the 
values of S given in table (IIVI) . where we also list the number of trials performed at each 
interface. 
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30 


1000 
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250 


1000 


18 


500 


1000 
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50 


1000 


10 


280 


1000 


19 


550 


1000 
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70 


1000 


11 


300 


1000 


20 


600 


1000 
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100 


1000 


12 


330 


1000 


21 


650 


1000 


4 


130 


1000 


13 


350 


1000 


22 


700 


1000 


5 


150 


1000 


14 


380 


1000 


23 


750 


1000 


6 


180 


1000 


15 


400 


1000 


24 


800 


1000 


7 


200 


1000 


16 


430 


1000 


25 


850 


1000 


8 


230 


1000 


17 


450 


1000 


26 


950 


1000 



Table IV: Interfaces and the number of trials per interface for the FFS sampling for the two dimen- 
sional Ising nucleation problem. 



The FFS calculation for the forward transition from A to B is straightforward. The 
flux <3>a through A from A is $^ = 1.5 x 10 -5 MC step -1 spin -1 and the forward rate 



constant k 



AB' 



2.8 ± 0.3 x 10 MC step spin : this is in good agreement with the 



value of 3.3 x 10 13 MC step 1 spin 1 computed for the same system by Sear et al. 



29fl. The 



computed forward rate constant does not depend on the reflecting wall position Sb'- This 



calculation also results in the function t + (S; A ), as described in section HVl In the reverse 
direction, we use the same interfaces and sample from X n to A as described in section [iVl 
We obtain the flux & B = 1.4 x 1CT 6 MC step -1 spin -1 and the backward rate constant 
k B ' A = 2.0 ±0.2 x 10~ 19 MC step -1 spin -1 . In this procedure, we also compute the function 
t_(S; A n ) j as described in section HV1 Combining the rate constants as in Eqs. lfTOl) and (fTTI ). 
we obtain pa = 7 x 10~ 7 and ps 1 = 0.999. By means of Eqs.(jlJ) and (JHJ) we finally obtain 
p(S) for 30 < S < 1000. Fitting this together with the distribution obtained by conventional 
sampling in state A (as described in section HvT) . we obtain the free-energy barrier. 




100 200 



S 

Figure 7: Free energy barrier for f3J = 0.65 and (3h = 0.05 calculated using FFS (dashed line; 
grey) and umbrella sampling (continuous line; black). Error bars are shown for both calculations. 



Figure [7] shows the results for the nucleation barrier, f3AG(S), in the range < S < 
300, computed as — In [p(S')]. The free-energy minimum at S ~ 20 indicates that for this 
supersaturation, the system has a small number of up spins even in the "down" state. The 
free-energy barriers, as obtained by umbrella sampling and FFS, are AG umbr = 24.5/cbT and 
AG FFS = 23k B T, respectively. These coincide within the error bars, which for both schemes 
are on the order of k B T. The computed barrier heights also agree remarkably well with the 
CNT prediction of 22k B T. 



IX. GENETIC SWITCH 



Our final test system is a biologically inspired non-equilibrium rare event problem: a 
model bistable genetic switch. This is a set of chemical reactions, representing protein- 



Figure 8: Schematic representation of our model switch, corresponding to Eq. [28l Two divergently- 
transcribed genes are under the control of a shared regulatory binding site on the DNA (the oper- 
ator). Each protein can bind, in homodimer form, to the operator and block the production of the 
other species. 



protein and protein-DNA interactions, as well as protein production and degradation, in 
a biological cell. The set of reactions shows two stable states, between which the system 
flips when simulated with stochastic dynamics. This is a particular case of the "exclusive" 
bistable genetic switch studied by Warren et al. [22|. The system does not obey detailed 
balance, and is therefore out of equilibrium. The set of chemical reactions which we simulate 
is given in scheme (|28| ). 



Reaction 


Rate 


Reaction 


Rate 




A + A ^ A 2 


k h k h 


B + B ^ B 2 


k { , k h 


(28a) 


+ A 2 ^ OA 2 


kon i k s 


+ B 2 ^ OB 2 


kon ; "'off 


(28b) 


-> O + A 


^prod 


-> O + B 


^prod 


(28c) 


OA 2 -> OA 2 + A 


^prod 


OB 2 -> OB 2 + B 


^prod 


(28d) 


A^0 


/' 


B -> 


f 


(28e) 



Our model switch is shown schematically in Fig. [H 

It consists of two genes, which encode proteins A and B. Proteins A and B can form 
homodimers A 2 and B 2 , as in Eq.(l28~k). The production rates for A and B depend on the 
state of the DNA sequence O, which is a regulatory site to which either A 2 or B 2 can bind. 
When O is free (not bound by either dimer), both genes can randomly be activated and 
produce either protein A or B with the same production rate /c pro d, as in Eq. (l28b ). When an 
A 2 dimer is bound to O (Eq. fl28b )). the production of B is blocked. Conversely, when a B 2 



dimer is bound to O (Eg. (125b)). the production of A is blocked. Both proteins can decay in 
the monomer form (accounting for active degradation processes and dilution in a growing 
cell), as in Eg.(l25e). We assume that transcription, translation and protein folding can be 
modeled as a single Poisson process, representing protein production. Clearly, when one 
species is abundant over the other one, many dimers of the majority species will be created, 
and the probability of finding one of them bound to O will be high. This effect will in turn 
lower the production rate of the minority species, leading to a stabilization of the state. If a 
rare fluctuation, however, is able to build up a substantial number of the minority species, 
these will in turn dimerize and bind to O, leading to a stochastic flip of the the switch. 



A mean field analysis carried out in [22] confirms this intuitive fact analytically: for 
suitable choices of the reaction rates, the system exhibits three fixed points: two symmetrical 
stable states, one rich in A and another rich in B, separated by one unstable state where 
the total number of A eguals the total number of B. The system can then be considered as 
a true bistable switch. 

We have chosen parameters such that the system is bistable and symmetric. Using the 
production rate k~ Iod as a time unit, and indicating by V the dimensionless volume of the 
system, we use: k{ = 5k pro aV, k^ = 5fc pro d (so that the eguilibrium dissociation constant 
for dimerization is Kjj = k^/kf = 1/V), k on = 5/c pro d, fc fr — ^prod (so that the eguilibrium 
dissociation constant for operator binding is = k s/k 0D = 1/(5V)), n = 0.3k. For 
simplicity, we will assume V = 1. The system is simulated with an event-driven Kinetic 



Monte Carlo algorithm 



32J which propagates the system according to the Chemical Master 



Eguation, thus accounting for the stochasticity arising from molecular discreteness and from 
the intrinsic randomness of reaction events. The simulation variables are the numbers of 
molecules n (copy numbers) of each chemical species. Briefly, in this algorithm, one selects 
at each simulation step a waiting time until the next reaction, and an identity for the next 
reaction, from the correct probability distributions. One then advances the simulation time 
by the chosen waiting time, executes the chosen reaction, and updates the copy numbers of 
the species involved in the reaction. 

A natural "order parameter" for the system is the difference between the total numbers 
of A and B proteins: q = X = + 2n^ 2 + 2n Q A 2 ~ ( n B + 2n B2 + 2n OB2 ). Since the system 
is symmetric, we know that $^4 = $b, k^B = &ba, an d therefore Pa = Pb = 0.5. As 
this system is out of eguilibrium, we do not sample a free-energy profile, but rather the 



non-equilibrium stationary probability distribution p(q) = p(A). 

To measure the switching rate and p(A), we run an FFS simulation with 12 interfaces, 
setting Ao = —27, X n = 27, and using N\ =10000 points at the first interface. The interfaces 
are positioned as shown in Table [VI 
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Table V: Interfaces and the number of trials per interface for the FFS simulations for the model 
genetic switch. 



We repeat the FFS sampling 10 times to obtain error bars. The result is Uab = ^ba = 
(8.66 ± 0.07) • 10 _6 fc~^ d . From the FFS calculations, we also obtain the function ^a(A) = 
iPa(q) as described in section IVTT1 and since the system is symmetric, we can obtain ^b(A) 
from iPa(X) by a simple inversion transformation. Combining i/ja{X) and V's(A), we arrive 
at p(A) = p(q) for —27 < q < 27, which is plotted in Fig. [9] (a scaling factor is applied to 
account for the different normalisation to the brute force results). The distribution is clearly 
bimodal and shows symmetric peaks whose positions correspond to the stable solutions of 
the mean field equations (|22|). As expected, a minimum in p(A) is observed for A = 
(unstable solution of the mean field equations). 

This system has a switching rate which is not exceedingly low, and we are also able to 
compute p(q) using a brute force simulation of length 2 ■ 10 9 k~^ od . The resulting stationary 
probability distribution is also shown in Fig. M Excellent agreement is obtained between 
the results of the FFS and brute force calculations. Because the system spends little time 
in the region between the two basins, this part of p(A) is hard to calculate accurately with 
the brute force run. The inset in Fig. [9] magnifies this region, showing the smooth profile 
produced by the FFS sampling. 




Figure 9: Probability distribution as a function of the order parameter A. The results are obtained 
both via long, steady state simulations (continuous line) and Forward Flux Sampling (circles). The 
region around |A| = can be accurately sampled only with FFS: the inset shows, on a logarithmic 
scale, a much smoother profile of the region close to the unstable steady state when FFS is used over 
Brute Force (BF). A scaling factor has been applied to the FFS results since they were originally 
normalised over —27 < A < 27 while the BF results are normalised over — oo < A < oo. 

X. DISCUSSION 

The key concept used here to obtain the stationary distribution in the unstable region 
between two stable states A and B is to add the contributions from the trajectories that start 
in A and go to B or return to A, and those that start in B and go to A or return to B (see 
Fig. [I]). These contributions can be obtained by performing one FFS calculation starting 
in state A and another starting in state B. For many rare event problems this is entirely 
possible - however, for systems where one state is very much more stable than the other, 
sampling the reverse transition (B — > A) may be computationally difficult, even with FFS. 
We have encountered this problem in the Ising nucleation example discussed here in section 
IVIIIi For equilibrium systems, this problem can be overcome by imposing an artifical stable 
state, as demonstrated here for the case of nucleation. However, this trick is not applicable 
for non-equilibrium systems. In general, in equilibrium systems the flux between any two 
state points is zero in steady state, while for non-equilibrium systems this need not be the 



case. In these non-equilibrium systems, the stationary distribution depends upon the full 
history of the trajectories. This, in general, prohibits the introduction of artifical boundaries. 
In particular, while for equilibrium systems detailed balance and microscopic reversibility 
dictate that the forward and backward transition paths have to occupy the same region in 
state space, for systems that are out of equilibrium the backward and forward trajectories 
do not have to coincide; indeed, in these systems cycles in state space can occur. We have 
recently demonstrated that the switching pathways of genetic switches can follow snch a 
scenario [Ufl. If the forward and backward transition paths form a cycle in state space, 
then it is conceivable that the artifical stable state "short cuts" the cycle and generates 
a wrong ensemble of points from which trajectories are initiated in the reverse direction. 
It may be possible to devise alternative techniques for sampling the reverse transition in 
non-equilibrium systems, and this will be the subject of future work. 

For the computation of free-energy barriers in equilibrium systems a wide range of nu- 
merical techniques is available The advantage of the scheme proposed here is that the 
free-energy can be directly obtained from an FFS simulation, obtaining simultaneously the 
rate constant, transition paths and free energy landscape. This is important because both 
the calculation of rate constants and the evaluation of free-energy barriers are computation- 
ally demanding, especially for large and complex systems. 

It has long been appreciated that free-energy barriers are critical quantities for under- 
standing rare events in equilibrium systems, such as nucleation and protein folding. However, 
the "barriers", or minima in the stationary probabilities, that separate steady states in non- 
equilibrium systems are equally important, because the rate of switching from one steady 
state to the next is proportional to the probability of being at the top of the "barrier" [221 ] . 
Some such "barriers" have recently been determined experimentally, including bimodal dis- 
tributions of protein concentrations for genetic switches like the one discussed in the previous 
34 l35| . To our knowledge, this technique is the first to be proposed for efficient 



section 



computation of stationary distributions for rare events in multi-dimensional non-equilibrium 
systems. This should prove useful for enhancing our understanding of a range of important 
non-equilibrium rare event processes, as well as improving the efficiency of computation of 
free energy landscapes in equilibrium systems. 
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Appendix 

In this appendix, we justify Eq. [6l Let us first imagine a very long brute force simulation 
trajectory which meanders around the basin of attraction of A, making occasional excursions 
towards B. We can divide each of these excursion into portions separated by successive 
crossings of interfaces Ao • . . A n . Consider the portion of an excursion between its leaving 
A and either reaching Ai or returning to A. We denote the distribution function (averaged 
over many excursions) for points visited during this portion fo(q) . Likewise, the distribution 
function (averaged over many excursions) for points visited after crossing Ai and before 
reaching either A 2 or A is denoted vi(q), and we can also obtain distribution functions Vi(q) 
for all interfaces < i < n. It is important to note that the Ui(q) are not normalised. In 
fact, the integral J dqv^q) contains information on the probability of an excursion reaching 
Aj. Since our entire ensemble of excursions can be divided up in this way, we can write the 
entire distribution function r + (g; A ) as the sum of contributions from all the portions of 
trajectories: 

ra-l 

r + (g;Ao) = 5>i(<2) ( 29 ) 

Now let us consider the FFS procedure. Let us imagine we have generated a collection 
of points at interface Aj. We fire Mj trial runs from this collection of points and continue 
each one until either A or A i+ i is reached. We plot a histogram 7i+(q; A«) of q values for the 
points in this ensemble of trial runs. We have proved before [12j that the distribution of 
these trial paths is identical to the distribution of corresponding portions of the "excursions" 
from A in a brute-force simulation, except that it is reweighted by a factor that depends on 



the probability of reaching A, from A - so that: 



Mo) 

P(A 4 |A C 



(30) 



We have also proved before 



12| that 



i-l 



P(A i |A ) = n p (Vil^ 
3=0 



(31) 



for i > (for % = 0, P(A |A ) = 1). 

Combining f l30l and ( 13T1) . we arrive at 

tt+(?; Ai) = 



nLo-p(\+ii^j 



(32) 



for i > and vr + (g; A ) = fo(?)- Rearranging Eq. [32] and summing over interfaces, we arrive 
at 

n— 1 n—1 i—1 

T+(g;A ) = J^i(g) = 7T + (g;A ) + ^7r + (g;Ai)JJP(A j+ i|A i ) (33) 

i=0 i=l j=0 

which corresponds to Eq. [6] 
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